Basic pipeline¶

Introduction¶

We'll get started importing scycle itself, and from there using single-cell RNA-seq from an Ewing Sarcoma cell line named CHLA9 [ref for original work]. scycle accepts AnnData objects as imputs. In this case, we internally download and the data from a loom file, which will also allow us to showcase the integration of the pipeline with scvelo and RNA velocity.

Pre-processing¶

Though samples can be pre-processed by other pipelines before using scycle, we implemented our own pre-processing pipeline. In the filter_cells stage, only quality control attempting to deal with low-quality cells is done. We also have integration with DoubletDetection to perform doublet removal inside of this step.

Then, the standard preprocessing recipe for scycle is called prep_pooling. This recipe includes an unconventional step that we find improves in the annotation of cell cycling: pooling (i.e. smoothing) cells using their nearest neighbors in a reference space. Otherwise, this step also performs library size normalization, log2-transformation and selects the most highly varying genes in the dataset.

An optional step in the pipeline that may improve downstream analysis, particularly finding the moment of cell division in the trajectory, is to normalize counts with reference to the mean total counts by partition. This is performed in a single-step in scycle in the pp module:

Dimensionality reduction¶

tl.dimensionality_reduction reduces the dimensionality of the data, by default, using Independent Component Analysis (ICA), as implemented in StabilizedICA. From the stable ICs, scycle then identifies the independent components related to cell cycle. Four possible dimensions are investigated: the two main known signatures of cell cycle, comprised of G1/S and G2/M, as well as 2 additional processes that are often but not always represented in the data. The Histones signature, which is comprised of histone genes, increases through S-phase and is stable through G2. The G2/M inhibitory (or G2-M-) signature related to a 'delayed' G2/M signal, which is still high at the start of G1 and whose transcripts start degrading later than transcripts forming the G2/M signature. This selection is performed by tl.dimensionality_reduction if method = 'ica' and find_cc_comps = True, which is default. It can be performed as a separate step invoking tl.find_cc_components.

In this example, we'll showcase an alternative dimensionality reduction methodology (method = 'self_consistent_CC), which starts from the approximation that uses the ICA followed by selection of cell-cycle components, and then uses the initial approximation of the cell cycle space to find genes related to cell cycle. These genes will be used as the reference for a principal component analysis, which will have a higher dimensionality space (n_comps-D) than the 3-4 D space found by the ICA. This has some benefits when calculating cell cycle trajectory later.

We can inspect our scycle results visually using the plotting module, but also by looking at the AnnData itself.

The previous steps have added multiple layers and variables to the AnnData. Here, we're going to focus on the obsm, which are dimensionality reductions. Feel free to take a look at other additions (e.g. chla9.uns['scycle'] will have the parameters for functions that have run)

In this case, X_dimRed is a 30-D PCA in the self-consistent cell-cycle genes, while X_pca_scycle is also a PCA with 3 PCs, which is used by the plotting module. scycle provides 2-D and 3-D visualizations using cc.pl.cell_cycle_pca and cc.pl.cell_cycle_pca3D

Trajectory and cell division¶

From this image, we can clearly see the trajectory, but also the effect of the partition counts normalization that segments the counts. We can get the trajectory from the principal elastic circle estimation using elpigraph, implemented in scycle as tl.trajectory:

Now we can look at the trajectory. We'll also plot it using the unormalized total_counts (total_counts_raw). We can still see a clear gradual increase in total counts, and the sudden drop-off at the moment of cell division (around node 6-7). Of note, when performing the partition normalization and using the self-consistent CC genes for dimensionality reduction, we often can find a "dimple" at the moment of cell division.

tl.cell_division will attempt to find the moment of cell division using total_counts (in this case, the normalized counts). It will set the direction of the trajectory as the direction of increasing counts, and the moment of cell division as the largest dip in total counts from the previous node. In this way, we can "break" the circle into a linear trajectory, and calculate pseudotime.

Pseudotime¶

We can see the distribution of cell over pseudotime, but now we can in fact plot many variables vs pseudotime in our data. For example, we can look at un-normalized counts:

Or the expression of some selected genes:

And, obviously, at the dynamics of cell cycle signatures over pseudotime. We could use pl.pseudotime_scatter again, but this is implemented into another call to pl.cell_cycle_scores.

tl.cell_cycle_phase uses the curvature of the trajectory to find inflection points in the expression of genes that define the cell cycle and set boundaries between cell cycle phases. It set the boundaries at the points in the trajectory where the curvature is the highest.

We can check how the cells were assigned into cell cycle phase 'categories' by looking at the pieplot:

or plotting using our friend pl.cell_cycle_pca setting col_var = 'cell_cycle_phase':

Cell cycle and RNA velocity¶

As scycle is built on top of AnnData, it is compatible with many other Python libraries that are built on the same platform, such as scanpy and scvelo. Since our example was loaded from a loom file, we have access to the unspliced and spliced matrices, and can look at RNA velocity in the cell cycle space

We can see in the stream plot a clear, orthogonal confirmation of inferred pseudotime directionality through RNA velocity.